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Research on transport, self-assembly and defect dynamics within confined, flowing liquid crystals 
requires versatile and computationally efficient mesoscopic algorithms to account for fluctuating 
nematohydrodynamic interactions. We present a multi-particle collision dynamics (MPCD) based 
algorithm to simulate liquid-crystal hydrodynamic and director fields in two and three dimensions. 
The nematic-MPCD method is shown to successfully reproduce the features of a nematic liquid crys¬ 
tal, including a nematic-isotropic phase transition with hysteresis in 3D, defect dynamics, isotropic 
Frank elastic coefficients, tumbling and shear alignment regimes and boundary condition dependent 
order parameter fields. 


I. INTRODUCTION 

As a state of soft condensed matter with interme¬ 
diate symmetries between highly ordered crystals and 
disordered fluids, nematic liquid crystals are both phe¬ 
nomenologically fascinating and commercially valuable. 
No longer are liquid crystals of interest only to those pro¬ 
ducing liquid crystal display technology; now scientists 
interested in microfabricated systems [T], microelectrome¬ 
chanical devices [2], composite materials [3], biosciences [4] 
and active gels [5] are exploiting the unique properties of 
liquid crystals in novel applications [6]. Interest in com¬ 
plex geometries (such as confining geometries nanocon- 
fined geometries [7], topological microfluidics [H |9] and 
colloidal intrusions [TOl [11]) require versatile mesoscopic 
algorithms that can account for non-trivial boundary 
conditions. Likewise research into “hypercomplex liquid 
crystals” m and self-assembly [T3l[T4] would benefit from 
efficient methods to simulate nematohydrodynamic baths 
for macromolecular and colloidal solutes. 

Such elaborate systems present a considerable chal¬ 
lenge for traditional particle-based numerical methods. 
Lattice Monte Carlo simulations have been very suc¬ 
cessful in simulating nematic liquid crystals [15] and con¬ 
tinue to be widely employed due to their computational 
frugality [ini fTT] . However, out-of-equilibrium dynam¬ 
ics and relaxation mechanisms require more computa¬ 
tionally costly methods. Off-lattice simulations of hard 
anisotropic particles and soft pair-potentials have played 
an important role in understanding generic liquid crys¬ 
talline phases UHlHo], but are limited to simple systems. 
Molecular dynamics simulations can account for molec¬ 
ular detail with a range of coarse-graining [2Ql m], in¬ 
cluding fully atomistic [22], generic molecules [23] and the 
mesoscopic approach of dissipative particle dynamics [2^ 
[25]. Even mesoscopic simulations can become compu¬ 
tationally expensive when large numbers of constituent 
particles are required and so are generally limited to sim¬ 
plified systems. 
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Investigating hypercomplex fluids or dynamics within 
demanding geometries calls for the continued develop¬ 
ment of versatile and computationally efficient coarse¬ 
grained algorithms. One mesoscopic simulation tech¬ 
nique that has shown promising capabilities in simulat¬ 
ing fluctuating hydrodynamics of isotropic solvents is the 
multi-particle collision dynamics (MPCD) algorithm 
[27] . MPCD has been used to simulate hydrodynamic 
interactions between macromolecules [28l [29] colloids [3Ql 
El], vesicles [32] and swimmers [33H35] . It has even 
been extended to simulate viscoelastic fluids [36] and 
electrohydrodynamics [37]. In this work, we propose an 
extension to the MPCD method to efficiently simulate 
fluctuating nematohydrodynamics {nematic-MPCD). 


II. METHOD 

Multi-particle collision dynamics algorithms forgo sim¬ 
ulating molecular-scale interactions between constituent 
molecules. Instead, the continuum description is discre- 
tised into many artificial, point-like MPCD particles that 
stochastically exchange momentum while respecting con¬ 
servation laws for mass, momentum and energy. This is 
sufficient to reproduce the hydrodynamic equations of 
motion on sufficiently long length and time scales. Meso¬ 
scopic MPCD algorithms can dramatically reduce com¬ 
putational costs compared to simulations that explicitly 
calculate molecular pair-potentials and are well suited to 
simulating flowing systems involving non-trivial bound¬ 
ary conditions [38l [39] . finite Reynolds numbers [40], and 
fluctuating hydrodynamics, which are ideal for moderate 
Peclet number systems [dT] [42]. 

Here, we develop a nematic-MPCD method to effi¬ 
ciently simulate fluctuating nematohydrodynamics, by 
assigning an orientation pseudo-vector to each MPCD 
point-particle and updating orientations through a local 
and stochastic nematic multi-particle orientation dynam¬ 
ics (MPOD) operator. Backflow and shear-alignment dy¬ 
namics are ensured by coupling the MPCD and MPOD 
operators. In § [ml we demonstrate that nematic-MPCD 
reproduces the necessary physical properties to simu- 
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late a nematic liquid crystal when the velocity and di¬ 
rector fields are coupled. In a very recent article, Lee 
and Mazza introduced an interesting hybrid, non-local 
MPCD method for liquid crystals [43]. The main differ¬ 
ence to our approach is that their particles carry a direc¬ 
tor field that is coupled to the fluid through a discreti¬ 
sation of the stress terms in a simplified Ericksen-Leslie 
formalism of nematohydrodynamics. 

In this section, we begin by reviewing a traditional 
Andersen-thermostatted MPCD algorithm that con¬ 
serves angular momentum. We go on to describe the 
implementation of the MPOD operator for nematic flu¬ 
ids and the two-way coupling between the director and 
velocity fields. Finally, we describe how potentially com¬ 
plex boundary conditions can be implemented. 


A. Traditional MPCD for Isotropic Fluids 


/cbT and is the average of the Nc random veloc¬ 

ity vectors in the cell during the instant of the colli¬ 
sion event. Randomly generating the from the equi¬ 
librium distribution /vei in the moving reference frame 
ensures that the algorithm is locally thermostatted [46] . 
The third term in the collision operator is a correction 
included to remove the angular momentum introduced 
by the collision operator 

Nc 

5L^ = E {^J ^ (% - i) } • (4) 
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Though the nematic-MPCD method does not strictly de¬ 
pend on this choice for coupling the velocity field to 
the director field is accomp lished by respecting this con¬ 
servation law (see § II C 2). 


The fundamental insight of MPCD algorithms is that 
continuous mass and momentum fields can be discretised 
into MPCD point-particles (labelled i). Each MPCD par¬ 
ticle possesses a position r^, mass rrii and velocity and 
which interact through multi-particle, near-equilibrium 
stochastic collision events within lattice-based cells (la¬ 
belled c) defined by a size a, population Nc^ centre of 
mass velocity centre ~ i-j )n moment of in¬ 
ertia ~ L/eL/c) point-particles 

in cell c relative to their centre of mass where 

- ry* ry* 

—i —i —cm,c* 

The MPCD algorithms consist of two steps. Each 
MPCD particle streams ballistically for a time 5t such 
that its position at time t ^ 5t becomes 

r. {t + 5t) = Ti (t) + Vi (t) 6t. (1) 

Multiple particles then undergo collision events, in which 
momentum is transferred between MPCD particles. To 
exchange momentum, the simulation domain is parti¬ 
tioned into cubic cells of thermally varying number den¬ 
sity pc = Nc/cl^ in d-dimensions. Discretising space 
into MPCD cells breaks Galilean invariance, though this 
can be remedied by randomly shifting the cell grid at 
each time step [44]. The collision operator is a non¬ 
physical exchange designed to be stochastic and also to 
conserve the net momentum within each cell c, 

Vi {t + St) = (t) + Si,c- (2) 

Many choices for the collision operator exist, which 
result in different versions of MPCD, including the origi¬ 
nal Stochastic Rotation Dynamics [26l \4E\ and a Langevin 
version of the algorithm [46]. In this work, we utilise the 
Andersen-thermostatted collision operator [46] 

where is a random velocity drawn from the Maxwell- 
Boltzmann distribution /vei (^, /cbT) for thermal energy 


B. Multi-Particle Orientation Dynamics for 
Nematic Fluids 

We now that propose nematic liquid crystals can be 
simulated via a nematic-MPCD algorithm by including 
an orientation field. 

Each MPCD particle is assigned an orientation u-^ 
while each cell acquires a tensor order parameter 

^ 

For a nematic fluid, the largest eigenvalue is the local 
scalar order parameter Sc of the cell and the local Frank 
director is parallel to the corresponding eigenvector. 

Orientations interact through a positive, globally spec¬ 
ified interaction constant U. In physical liquid crystals, 
the energy U represents inter-molecular interactions and 
will be a non-constant function of temperature or molec¬ 
ular details such as nematogen dimensions and density. 
In this nematic-MPCD algorithm, the interaction con¬ 
stant U is the simulation specified energy that governs 
the local evolution of orientations. Taking inspiration 
from the Andersen-thermostatted MPCD collision oper¬ 
ator, we implement a stochastic multi-particle orientation 
dynamics operator for orientation. The essential require¬ 
ments are that the MPOD operator must be local and 
near equilibrium^ with no gradient terms in the collision 
operator. Therefore, we propose the orientation collision 
event 

Ui {t + St) = 4'c (u, (i)) , (6) 

where the multi-particle orientation operator Tc gen¬ 
erates a random orientation u- {t + 5t) drawn from 

the equilibrium probability distribution /ori Q (t)^ 

about the local director (t) calculated from the tensor 
order parameter. 




3 


1. Maier-Saupe Distribution: 


1. Shear Alignment: Veloeity^ Orientation Coupling: 


As in the traditional Andersen-thermostatted MPCD 
algorithm, the multi-particle orientation operator de¬ 
pends on the condition of local, near-equilibrium statis¬ 
tics. In this work, we assume that the local equilibrium 
distribution for the orientation field obeys the Maier- 
Saupe self-consistent mean-field theory and so is an ex¬ 
ponential function of Un = Ui ’ He' 


The nematogens respond to the flow field’s vor- 


ticity u 


Vv — (Vv) ^ /2 and shear rate D 


+ {^31) /2 (Fig, 


8 insets) by obeying the discre- 


tised Jeffery’s equation for a slender rod 


St 


Xhi ;£)] , (9) 


/ori {u, S„ri,) = (7) 

where ^ is a normalisation constant, P = 1/k^T and each 
cell’s mean-field interaetion potential is 

wmf,c = —UScU^ + — {Sc — 1 ). ( 8 ) 


where A is the bare tumbling parameter and xhi is the 
shear eoupling eoeffieient^ a simulation parameter that 
tunes the alignment relaxation time relative to St. 

For the rotation of an individual ellipsoidal particle 
subject to shear flow, xhi = 1- When the shear coupling 
coefficient is set to zero (xhi = 0 ) there is no coupling of 
the director to the velocity field. 


The second term does not depend on and so the distri¬ 
bution of Un is determined by . When the scaled 

energy /SUSc is small, all orientations are equally likely 
but when PUSc is large the distribution becomes sharply 
oriented about n^. 


2. Generating the Maier-Saupe Distribution: 


When /3USc ~ 1, a Metropolis algorithm for wmf,c 
generates the random orientations. However, the dis¬ 
tribution can be more efficiently approximated in the 
limits of pUSc ^ 1 and /SUSc 1. In the strong mean 
field limit ^USc ^ 1 , fori is sharply centred about 
such that u‘^ = cos^ On ~ 1 — which means that the 
distribution for On can be approximated as Gaussian 
/ori ^ The Gaussian approximation is used 

when jSUSc > 5. On the other hand, when jSUSc 1 
the exponent can be expanded and the cumulative 
distribution function of /ori can be approximated as 
W = Un ^ pUScul/Z. Random 

values of Un can be generated through the transfor¬ 
mation Un = njPUSc — where K{r) = 

r 11 / 2 \ 

3r(/3[/5c)^+ ) and 

r G [0,1]. This expansion is used when pUSc < 0.5. 


C. Two-way Coupling 


Goupling between the director and fluid flow is cru¬ 
cial for reproducing nematoh ydrodyn amics since flows 
can rotate the nematogens (§ |II C 1[ ) and the rotation 
of nematogens in turn pro duces h ydrodynamic motion, 
referred to as baekflow (§ IIIC2D . We model the cou¬ 
pling to be superimposable overdamped torques. 


2. Baekflow: Orientation-^ Veloeity Coupling: 

The nematogens are implicitly envisioned as rotating 
through the viscous fluid that they themselves represent, 
and hence they experience viscous torques. Assuming 
that the rotational motion is overdamped, we model this 
as a rotational Stokes drag. In the nematic-MPGD algo¬ 
rithm, backflow coupling is accounted for by balancing 
this drag on the nematogens via transferring an equal 
and opposite change in angular momentum to the veloc¬ 
ity collision operator. 

The magnitude of the rotational Stokes drag felt by 
the nematogens is T- = x u -. The nematic collision 
causes the MPCD point-particles to change their orien¬ 
tation by JMcoI i — Hi {t S- St) — {t) over the time step 

Jt, and so results in a collisional contribution to the drag 
torque 

Ecol.i = X (10) 

Likewise, the torque on the fluid due to the shear align¬ 
ment is Enpi = iRHi X Su^i JSt from Eq. 1^ If any 
additional torques due to external fields wer^o be in¬ 
cluded they too would need to be included in the net 
torque. 

To balance these torques with the hydrodynamic drag, 
the opposite of the net change in the angular momentum 
— Yin"" is transferred to the linear 
momentum portion of the algorithm. The MPCD colli¬ 
sion operator (Eq.|§ is thus modified to account for 
liquid crystal backflow becoming 

= £. - + (f ^ - «J) X (11) 

In this way, the total angular momentum of the system 
is conserved and the orientation-velocity coupling is ac¬ 
counted for. By setting 7 ^ = 0, the transferred angular 
momentum of each particle is zero SC^ = 0 and this cou¬ 
pling can be turned off. 
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FIG. 1: Nematic-isotropic phase transition. Simulation 
parameters from §211 are used with no shear coupling, 
Xhi = 0. Simulations in 3D exhibit discontinuous 
isotropic to nematic transitions, regardless of whether 
the local Sc (solid lines) or global S (dashed lines) order 
parameter is used. The transition is second order in 2D 
when the local Sc is used but becomes discontinuous if 
S is used. The nematic-isotropic transition agrees 
qualitatively with the Maier-Saupe self-consistent 
mean-field theory (MS; dotted lines). Inset shows a 
typical snapshot of the isotropic disordered state (left) 
and nematic ordered state (right). 


or a single preferred direction can be specified. If no pre¬ 
ferred direction is specified then the boundary is said to 
be non-anchoring. 


E. Units and Chosen Simulation Parameters 

Values are expressed in MPCD simulation units — 
time, mass, energy and length are given respectively by 
time step 6t, par ticle ma ss m, thermal energy /cbT and 
cell size a = jm. The new MPOD parameters 

are also stated these units. The interaction constant U 
has units k^T^ while the rotational friction coefficient 
has units k^T 6t. Both the bare tumbling parameter A 
and the shear coupling coefficient xhi are dimensionless. 

Except when otherwise stated, the simulations pre¬ 
sented in this manuscript vary input parameters about 
the following set of values: In this manuscript, simula¬ 
tions are carried out in 2D (d = 2) for a system of size 
V = 50^ with periodic boundary conditions and a mean 
number density of p = (Nc) = 20. The MPCD parti¬ 
cles are randomly initiated with positions from a uniform 
distribution, velocities from the Maxwell-Boltzmann dis¬ 
tribution and aligned nematic orientations. Parameter 
values are chosen to be m = 1 , /cbT = 1 , dt = 1, a = 1, 
U = 15, Jr = 0.01, A = 2 and xhi = 1- 


III. RESULTS 


D. Boundary Conditions 


One of the advantages of particle-based hydrodynamics 
solvers is that complex and mobile boundary conditions 
can be implemented. For this reason, the nematic-MPCD 
may be well-suited to nematic fluids confined within mi¬ 
crofluidic devices [H [9] and to simulating colloidal-liquid 
crystals pTl fTT] and hypercomplex liquid crystals [12]. 

The effect of boundaries on positions and velocities are 
implemented in the standard manner. Periodic bound¬ 
ary conditions are implemented by wrapping the MPCD 
particle positions. Lees-Edwards boundary conditions 
are used to introduce simple shear flows across peri¬ 
odic domains [47]. No-slip walls are simulated by imple¬ 
menting bounce-back boundary conditions with phantom 
particles [48H5Q] . 


The boundary conditions also set the easy direction 
describing the preferred orientation of the liquid crys¬ 
tal director at a surface. During a bounce-back collision 
event with the surface, the orientation of the imping¬ 
ing nematic-MPCD particle is set parallel to the surface’s 
easy direction. This anchoring is not strong, as will be 
seen in §211 ^ For homeotropic boundary conditions, 
the easy direction is normal to the surface. For planar 
boundary conditions, the easy axis is parallel to the sur¬ 
face. In this case, all in-plane directions can be equivalent 


Having described the implementation of the nematic- 
MPCD algorithm, we now characterise the resulting 
properties of the liquid crystal. We first consider how 
the isotropic-to-nematic phase transition depends on the 
simulation parameters, particularly the heuristic shear 
coupling coefficient and number density of MPCD par¬ 
ticles. We measure the nematic-isotropic hysteresis and 
explore the dynamics of the defect annihilation rate as 
the system orders. Elastic free energy drives defect anni¬ 
hilation and we measure the isotropic Frank elastic coef¬ 
ficients to be a linear function of the interaction constant. 
The response of the isotropic phase to an ordering wall 
is characterised. 


A. Nematic-isotropic transition 

When /3I/ is small, the nematic-MPCD algorithm ex¬ 
ists as an isotropic fluid state with a small global order 
parameter S (Fig. inset-left). When pU is large a ne¬ 
matic state is formed (Fig. inset-right). Maier-Saupe 
self-consistent theory predicts that the nematic-isotropic 
transition is first order (Fig. [^. Although the nematic- 
MPCD algorithm assumes near-equilibrium and so uses 
the Maier-Saupe distribution on the local cell level, the 
scalar order parameter and directors are spatially varying 
fields rather than mean-field values. 
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FIG. 2 : Global order parameter as a function of simulation parameters varied about the simulation values from 
§ [irEl (red squares). The top row shows simulations in the absence of shear coupling, xhi = 0. The bottom row 

shows simulations with full coupling, xhi = 1 - 


In 3D systems of large enough size, periodic boundary 
conditions and no shear coupling, the nematic-MPCD 
algorithm does exhibit a strongly first order nematic- 
isotropic phase transition (Fig. Q. In these simulations, 
the nematic fluid is initialised in the nematic state and 
resides in a periodic cube of size 50^. The system dis- 
continuously jumps from zero to a global scalar order 
parameter S'* = 0.860 ± 0.003 at = 4.20 ± 0.05. 

In 2D the nematic-isotropic transition is expected 
to become a Kosterlitz-Thouless-type transition [5T1 [52] . 
The present simulations demonstrate that the transi¬ 
tion is no longer first order, increasing from zero at 
[ydl/]* = 4.1 ± 0.1 in a 50^ system (Fig. [^. The sec¬ 
ond order nature of the nematic-isotropic transition is a 
direct result of the nematic-MPGD’s ability to accom¬ 
modate spatialtemporal varying fields. Future studies 
should more fully characterise the nature of the nematic- 
isotropic transition in 2D. 


1. Global vs. Local Scalar Order Parameter: 


By replacing the local scalar order parameter Sc of each 
cell c with the system’s globally determined order param¬ 
eter S in each cell’s local mean-field interaction poten¬ 
tial wmf c (Eq. [^, the 2D transition becomes first order 
(Fig. §. The order parameter curve remains relatively 
unchanged except near the phase transition. The tran¬ 
sition from the isotropically disordered state is retarded 
compared to the spatially varying case that uses the local 
order parameters but suddenly jumps to S'* =0.51±0.01 
at = 5.00 ±0.05. 


2. Variation of Simulation Parameters: 

In order to assess the impact of varying simulation 
parameters on the nematic-isotropic transition, we ini¬ 
tially omit the velocity^orientation coupling by setting 
Xhi = 0 in Eq. (Fig. |2] top row). We consider varying 
time step 5t^ mass m, temperature rotational fric¬ 

tion coefficient 7 ^^, bare tumbling parameter A and mean 
number density p. In the zero-coupling limit. Fig. 
(top row) shows that none of the MPCD simulation pa¬ 
rameters have a significant effect on the nematic order¬ 
ing. Only the mean number density has an observable 
affect on the curve. At extremely low mean number den¬ 
sities, the transition occurs at a slightly larger interaction 
constant. It should be noted that when an individual 
nematic-MPCD particle is alone in an MPCD cell nei¬ 
ther its velocity nor its orientation are altered. 


3. Impact of Coupling Fluctuating Hydrodynamics: 

When the shear coupling coefficient xhi is zero the 
global scalar order parameter S rises from zero in the 
isotropic phase to S' = 1 in the PU ^ 00 ordered limit. 
This is no longer true when xhi 7 ^ 0 (Fig. bottom 
row). As the hydrodynamic coupling is restored by in¬ 
creasing Xhi, the value of the scalar order parameter de¬ 
creases for a given interaction constant. This occurs be¬ 
cause fluctuations in the velocity field introduce an ad¬ 
ditional source of noise through Eq. E when Xhi 7 ^ 0. 
These fluctuations reduce the order in the director field 
and move the system away from the fully ordered state 
of S = 1. With full coupling, only the rotational friction 
coefficient jr is seen to have no impact on the S curve 
(Fig.i bottom row). This is because jr controls the ro¬ 
tational relaxation dynamics and does not influence the 
equilibrium state. 
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FIG. 3: Nematic-isotropic phase transition as a function 
of average number density for 1/ = 15. Simulation 
parameters from § EH are used. Without shear 
coupling (xhi = 0), the global order parameter S is 
relatively constant but when xhi = 1 the system 
exhibits a density dependent second order phase 
transition. 


When Xhi = 0, the mean number density is the only 
simulation parameter seen to have any observable effect 
on the isotropic to nematic transition and then only at 
extremely low values (Fig. top row). At the lowest 
number density the transition is less sharp and occurs 
at a slightly higher interaction constant PU. While S 
depends weakly on p when xhi = 0 (Fig. [^, it is a 
strong function of number density when xhi = 1- Fig. 
shows the global scalar order parameter as a function of 
density for U = 15. When the shear coupling parameter 
Xhi is set to zero, the system remains in the nematic state 
even at quite low densities. On the other hand, when 
coupling is included, S increases from zero with mean 
number density. In fact, the order parameter in Fig. 
exhibits a continuous transition and the nematic-MPCD 
algorithm possesses a nematic-isotropic transition as a 
function of density when xhi = 1- 

Since there is a nematic-isotropic transition as a func¬ 
tion of density (Fig. [^ , it is clear that the shear coupling 
coefficient has a larger effect at lower number densities 
than it does at larger densities. Fig. shows the strong 
interaction limit of S (measured at /3U = 100 and 500) 
for various densities as a function of coupling. For a low 
mean number density of p = 5, Fig. shows that the 
strong limit drops from limpt/^oo S ^ 1 when xhi = 0 to 
only 0.038 ± 0.003 when the algorithm is fully coupled. 
Fluctuations are pronounced because of the small number 
fluctuations of particles in each MPCD cell. By increas¬ 
ing the mean number density p, the continuum limit is 
approached and fluctuations become less severe. When 
p = 20 and the algorithm is fully coupled (xhi = jO ^ 
the strong interaction limit is S' = 0.80 ± 0.01 (Fig. Q. 


FIG. 4: Increased shear coupling (xhi) reduces the 
scalar order parameter. The order parameter is 
measured in the highly nematic phase at pU = 100 and 
500. The other simulation parameters are given in 


§ HE 


Throughout this work, we set the mean number density 
p = 20, though a lower density may suffice in many situ¬ 
ations. 

When the algorithm is fully coupled with xhi = 1, the 
tumbling parameter can also increase the susceptibility 
of the order field to velocity fluctuations through Eq. 
This can be seen in Fig. (bottom row). When A = 5, 
fluctuations in the shear rate D fully disorder the system. 


4- Hysteresis: 


Hysteresis is expected in 3D due to the first order na¬ 
ture of the nematic-isotropic transition. By comparing 
3D nematic-MPGD simulations initialised with the direc¬ 


tor field in the isotropic state (as in in § III A ) to those 
initialised in the nematic state, a striking hysteresis loop 
is observed in Fig. The interaction constant, [311^ is 
fixed throughout the duration of individual simulations. 
At these system sizes, the width of the hysteresis is mea¬ 
sured from Fig. to be /3AI/* = 0.70 ± 0.03 and the 
difference in order parameters at the transition points is 
A5'* = 0.20 ±0.04. 


B. Defect Annihilation Dynamics 

The process of transitioning from the isotropic to ne¬ 
matic phase discussed in § ME is controlled by the 
dynamics of topological defects. Though the increas¬ 
ing interaction constant U generates local order along 
a spontaneous direction neighbouring regions may 
break symmetry along any other direction. Therefore, 
many ±1/2 topological defects rapidly emerge from the 
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FIG. 5: Hysteresis loop in the isotropic/nematic 
transition. Simulations are initialized in the isotropic 
(circles) or nematic state (triangles). Simulation 
parameters from § HE are used in 3D. 


disordered director field. Pairs of oppositely charged de¬ 
fects must approach each other and annihilate for global 
ordering. 

Since the 2D number density po = 0.0080 ± 0.0005 
of defects is initially quite high, the annihilation rate 
i?D = —Pd is large but falls rapidly (Fig. inset a). 
As the density decreases, the average separation between 
topological defects increases and annihilation events be¬ 
come less frequent (compare Fig. |6a| showing an example 
system at t = 40 to Fig. |6b| showing the same system at 
t = 400). A variety of scaling relations for the annihila¬ 
tion rate Rd (t) ^ have been put forward. Mean- 

field arguments predict z/ = 1, purely diffusive kinetics 
suggest u = 0.5 and scaling arguments give u = 6/7[53]. 
Furthermore, the scaling law possesses short-time log¬ 
arithmic corrections [5^ [55] . When measured on short 
times between t G [lO, 10^], the nematic-MPCD annihi¬ 
lation rate appears to decay as z/ = 0.74 ± 0.02 (Fig. 
but the exponent increases to z^ = 0.83 ± 0.04 when eval¬ 
uated over t G [lO^, 10^] (Fig. [^, which is in agreement 
with the u = djl scaling prediction. 


C. Prank Elastic Coefficients 

We have considered how the nematic state arises from 
the isotropic state. Let us now consider the nematic re¬ 
sponse to distortions in the director field. Gradients in 
the director field n lead to the free energy density per unit 
volume / = Kspiay (V ■ nf /2 + iiltwist (zi ■ V x nf /2 + 
A^bend (zt X (V X n))^ /2. Splay, bend and twist deforma¬ 
tions are illustrate in Fig. [g insets. Since distortion 
are typically large compared to molecular length scales, 
the Frank elastic coefficients ATspiay, A'twist and A'bend are 



(a) Director field at t = 40St. (b) Director field at 

t = 400(5t. 



FIG. 6: Number density of topological defects as a 
function of time with and without shear coupling. 
Inset (a): Annihilation rate of defects Rjy. Simulation 
parameters from § IIE are used with a 2D system size 
of 500 X 500. In the defect maps, colour denotes local 
scalar order parameter and defects are mapped with red 
circles marking —1/2 defects and blue pentagons 
marking +1/2 defects. 


macroscopic material properties. 

One technique for obtaining the Frank coefficients from 
particle-based simulations is to measure the equilibrium, 
orientational fluctuation spectrum [561159] . In reciprocal 
space, the tensor order parameter for each wave vector is 
£ it) = P~^ + {dUiUi - i) exp {ik ■ rj. We work 

in a varying director-based coordinate system, in which 
= [0,0,1] and the wave vector is in the 13-plane, i.e. 
k = [ki^O^ks]. In this coordinate system, the equipar- 
tition theorem [56] relates the the low \k\ limit of the 
orientational fluctuations to the Frank coefficients 


Qas {k) Qas {-k)) = - 


9 SVkBT 


4 + A^bend^s 


( 12 ) 


for a = 1,2 (splay, twist). Through Eq. 12, the Frank 


coefficients may be determined as fitting parameters of 
the fluctuation spectrum in reciprocal space. A large 
system size of V = 30 x 30 x 30 and N = 5.4 x 10^ 
MPGD particles are used in the following simulations to 
ensure sufficient statistics for many near-zero \k\ values 
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FIG. 7: Frank elastic coefficients for splay, twist and 
bend as a function of interaction constant U. Simulation 


parameters from § |II E| are used for a system of size 
30 X 30 X 30 {N = 5.4 x 10^ nematogens) with xhi = 0. 
Insets a,b,c respectively depict splay, twist and bend. 


and accurate fits. 

The resulting Frank coefficients of the nematic-MPCD 
fluid are shown in Fig. [T] Although splay, twist and 
bend deformations may possess differing coefficients in 
some physical systems, simple scaling suggests that all 
three elastic constants are of order ^ U/a and theo¬ 
retical considerations of the Maier-Saupe self-consistent 
model[60] predict Ki = pUS‘^ [1 ^ Ci] / 6 ^ where i = 
{splay, twist, bend} and is a characteristic interaction 
distance. The different constants Ci depend on the 
molecular details and higher moments of the orienta¬ 
tion distribution [60]. The nematic-MPCD simulations 
ostensibly exhibit isotropic elasticity. This is expected 
because, in the limit that the rod length is small com¬ 
pared to the interaction length, the constants Ci are 
safely neglected and the Frank coefficients are predicted 
to converge [60]. Since the nematic-MPCD algorithm sim¬ 
ulates point-like nematogens with a characteristic inter¬ 
action length equal to the finite cell size, the one-constant 
approximation applies. 

In agreement with simple scaling and the Maier-Saupe 
self-consistent predictions, the measured elastic coeffi¬ 
cients for the nematic-MPCD algorithm are linear with 
respect to the interaction constant U (Fig.[^. Together, 
they are fit to = (113 ± 5) for the parameters in 
Fig.ig 


D. Tumbling and Shear Alignment 

Thus far, we have considered quiescent nematic flu¬ 
ids. We now turn our attention to flowing systems. Mi¬ 
croscopically, the director field is influenced by shearing 
flows through Eq. and as described schematically in 


Fig.m inset. 

In the infinitely dilute limit of a suspension of 
spheroidal particles, the bare tumbling parameter is a 
geometrical entity that can be cleanly related to the par¬ 
ticle aspect ratio p by A = (p^ — l) / (p^ + l), which 
goes to unity as p ^ oo and is zero for spheres (p = 1 ). 
However, interactions between nematogens in a nematic 
fluid allow the actual tumbling parameter to deviate 
from the isolated-slender-rod value and distributions of 
molecules can exhibit effective tumbling parameters that 
are larger than unity. Such fluids are referred to as 
aligning-nematics because there is a stable alignment an¬ 
gle, the Leslie angle Ol, between the director and shear 
field. 

By considering a Fokker-Planck equation for the prob¬ 
ability distribution of orientations. Archer and Larson [61] 
found that the flow tumbling behaviour of ellipsoidal par¬ 
ticles with A = (p^ — l) / (p^ + l) is determined by the 
tumbling parameter 


A' = A 


155' + 4854 + 42 
1055 


(13) 


In the nematic-MPCD algorithm, A is the specified sim¬ 
ulation parameter for the bare tumbling parameter, the 
magnitude of which can be s et la rger than unity. We shall 
see that A' as given by Eq. ^ is the resulting tumbling 
parameter of the nematic-MPCD algorithm. In Eq. |13[ 
54 is the fourth moment of the Maier-Saupe probability 
distribution. The distribution can be written as an ex¬ 
pansion of orthogonal Gegenbauer polynomials (x) 
in d-dimensions as 


^ 4.^+1 (d- 2 ^ 

fori (U, S,ri) = J2 (14) 

where the moments are 821 = /• In 3D, the 

polynomials are Legendre polynomials, while they are 
Chebyshev polynomials in 2 D. The first even moment 
is the scalar order parameter S = 82 = {u^ ~ 2 ) 

representing the variance of the alignments about the 
director, while 54 = 84 is the next non-zero moment. 


1 . Tumbling Nematic: 

When A' < 1, the nematogens continuously revolve or 
tumble. The tumbling period is set by the Jeffery orbits 
to be 


2tt 

Xhi7\/1 - A'2 ’ 


(15) 


where 7 is the shear rate. 

Using Lees-Edwards boundary conditions m to es¬ 
tablish a shear rate 7 = 0.01 across a periodic channel 
of height L = 50, we measure the tumbling period as 
a function of tumbling parameter A' (Fig. [^. The pe¬ 
riod is relatively small when A' is small and varies very 
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Tumbling Parameter, A' = A (155 + 48^4 + 42) /1055 Tumbling Parameter, A' = A (155 + 4854 + 42) /1055 


FIG. 8 : Jeffery periods corrected for the tumbling 
parameter and non-unity shear coupling coefficient. 
Simulation parameters from § ini are used with 
U = 20 and orientations initialized in the nematic state. 
The insets show the rotational and extensional 
components of Eq. 


FIG. 9: Leslie angles Ol corrected for the tumbling 
parameter. Simulation parameters from § im are used 
with U = 20 and orientations initialized in the nematic 
state. Ericksen diagram is inset and shows stable and 
unstable orientations for Ol. 


little as a function of tumbling parameter. However, as 
the tumbling parameter increases, the period increases 
rapidly and diverges as A' ^ 1. The simulated tumbling 
periods are found to be in good agreement with Eq. |15| 

The tumbling period does not depend on the rotational 
friction coefficient (Fig. [^. This is expected both 
from inspection of Eq. and from the realisation that 
the differential drag by the shearing ffow is what rotates 
the rod. 


2. Shear-Aligning Nematic: 



Interaction Constant, /3U 


When the magnitude of the bare tumbling parameter 
A is set so that |A'| is larger than unity, the nematogens 
do not tumble but rather align with the shear. For these 
tumbling parameters, Eq. has the solution 


Good agreement is found between Eq. |16| and the sim¬ 
ulations using Lees-Edwards boun dary conditions when 
the tumbling parameter (Eq. |13[ ) is greater than unity. 
As the tumbling parameter tends to 1 +, the Leslie an¬ 
gle approaches zero. In this limit, the nematogens orient 
along the ffow direction. When A' ^ 1, the Leslie angle 
of the nematic-MPCD ffuid approaches 7 r /4 as predicted 
by Eq. [T^ 


FIG. 10: Wall-induced ordering. Inset shows 
exponential decay of order from homeotropic wall 
towards the isotropic bulk, as characterised by the 
coherence length C Coherence length diverges as the 
nematic-isotropic transition is approached. Simulation 
parameters from § 


HE 


are used with 7 ^ = 1 . 


E. Wall-Induced Ordering 

Confining walls affect the nematic ordering. In the 
isotropic state, anchoring can cause ordering in the vicin¬ 
ity of the walls. We consider a 2 D nematic-MPCD 
ffuid confined between two no-slip plates separated by 
L = 100. The plate at ^ = 0 enacts homeotropic bound¬ 
ary conditions, which order the nematic ffuid. The plate 
at ^ = 1 / is a non-anchoring boundary, which does not 
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set a condition for u^. 

When the interaction constant is much l ess tha n the 
nematic-isotropic transition value (§ III A), the 

order decreases to the isotropic state far from the wall. 
As the interaction constant is increased, the value of the 
scalar order pa ram eter Sq{U) = S' ([/, ^ = 0) at the wall 
increases (Fig.p^ inset). This signifies that the anchor¬ 
ing is not infinitely strong and is strongly effected by the 
value of U. 

Additionally, the order extends further into the bulk 
fluid as U increases. The characteristic distance the 
order ext ends from the wall is a coherence length ^ 
(Fig. ra. One can predict that the order decays as 
S (1/, y) = So {U) by considering the total free 

energy functional to be the highest order term in the 
Landua-De Gennes free energy and the deformation free 
energy. The coherence length is a function of the elas¬ 
tic constant and the distance from the transition, ^ (x 

^ 3Kspi^y+^Ktwist ^ nematic-isotropic transi¬ 

tion is approached, the coherence length diverges. The 
nematic-MPCD simulations accurately reproduce the ex- 
ponential decay far below the transition point (Fig. |1Q[ 
inset). 

It was seen in § III C that Ki U so the coherence 
length takes the form 




cipU 

1 - U/W 


1/2 


(17) 


The theory captures the rapid growth of the coherence 
length near the nematic-isotropic transition but go es t o 
zero as 1/ ^0, while the simulations do not (Fig.[^. 
The coherence length of the nematic-MPCD does not 
go to zero because MPCD algorithms are not able to 
resolve material properties on length scales compara¬ 
ble to the cell size a. If a second fitting parameter 
C 2 = (1.505 ± 0.005) a is included as in Fig. 


10 


to ac¬ 


count for this discretisation effect, then Eq. Ill well- 
represents the divergence of the coherence length in the 
isotropic phase. 

In the nematic phase, the order still decreases expo¬ 
nential from So but decays to a non-zero bulk value 
(Fig. |10[ inset). Except near the nematic-isotropic tran¬ 
sition, the order parameter falls steeply to its bulk value 
over a length scale comparable to a single MPCD cell. 


MPCD uses traditional Andersen-thermostatted MPCD 
with conservation of angular momentum to integrate the 
velocity field and a novel multi-particle orientation dy¬ 
namics (MPOD) collision operator to progress the direc¬ 
tor field. By stochastically drawing orientations from the 
local Maier-Saupe equilibrium distribution, the MPOD 
operator updates the orientations without numerically 
evaluating gradients. In addition, the two-way coupling 
between the MPCD and MPOD operators represents 
backflow and shear-alignment. We have shown that this 
nematic-MPCD algorithm reproduces the essential phys¬ 
ical properties of a simple nematic fluid, such as the 
nematic-isotropic phase transition, topological defects, 
Frank elasticity and shear alignment. 

The nematic-MPCD algorithm holds much promise 
as a tool for simulating nematohydrodynamics, but fu¬ 
ture studies should carefully investigate the anchoring 
strength (since modifications to the no-slip conditions 
were required in traditional MPCD [48l[5Q] i and work to¬ 
wards kinetic theories to quantitatively predict the ne¬ 
matic material properties as a function of simulation pa¬ 
rameters. Though simple, the algorithm holds exciting 
potential for simulating a wide variety of soft matter 
systems. For example defect dynamics within topolog¬ 
ical microfluidic devices [9] or porous media [62] could be 
modelled, exploiting the ease with which the algorithm 
can handle complicated confining geometries. It would 
also be of interest to consider dispersed nanoparticles, 
carbon fibres [63] or swimmers [64] within a liquid crystal 
host, and it is relatively easy to imagine that general¬ 
ized Maier-Saupe theories [65] could be implemented to 
in the MPOD collision operator to simulate cholesteric 
or biaxial liquid crystals. 
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